#=====================================================================
# 2020/07/27
# Geothermal prediction dataviz
# Fumiya Uchikoshi, uchikoshi@princeton.edu
#=====================================================================
rm(list=ls())
######################################################################
# Loading packages
######################################################################
library(tidyverse) 
library(memisc)
library(dplyr)
library(tidyr) 
library(xtable)
library(ggplot2)
library(ggthemes)
library(magrittr) #set_colnames
library(rlang) #eval_bare

######################################################################
# Figure 1
######################################################################
df <- read.csv("../Data/prevex.csv",fileEncoding = "Shift_JIS",header=TRUE)%>% 
  mutate(prevexpyeardum=case_when(
    prevexpyeardum == 0 ~ "Not explored (n=582)",
    prevexpyeardum == 1 ~ "Explored (n=67)"
  ),
        potmax=exp(logpot))
    
table(df$prevexpyeardum)

ggplot(df, aes(x=prevexpyeardum, y=potmax,color=prevexpyeardum)) +
  geom_violin(trim=T,fill="#cfcfcf",linetype="blank")+xlab("")+
  ylab("Geothermal Activity Index")+theme_few()+theme(legend.position="none")+
  geom_hline(yintercept=5.86,linetype="dashed",colour="black")+
  annotate("text",label="Cutoff (=5.86)",x=0.65,y=15)+
  stat_summary(geom="pointrange",fun.y = mean, fun.ymin = function(x) mean(x)-1.96*sd(x), fun.ymax = function(x) mean(x)+1.96*sd(x), size=1,alpha=.5)
ggsave(height=4.5,width=6,dpi=200, filename="../Results/Figure1.pdf",  family = "Helvetica")
